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Abstract. We propose a method for the stabilisation of a stack of parallel vortex 
rings in a Bose-Einstein condensate. The method makes use of a "hollow" laser beam 
containing an optical vortex. Using realistic experimental parameters we demonstrate 
numerically that our method can stabilise up to 9 vortex rings. Furthermore we point 
out that the condensate flow through the tunnel formed by the core of the optical vortex 
can be made supersonic by inserting a laser-generated hump potential. We show that 
long-living immobile condensate solitons generated in the tunnel exhibit sonic horizons. 
Finally, we discuss prospects of using these solitons for analogue gravity experiments. 

PACS numbers: 03.75.Lm, 03.75.Kk, 04.70.Dy 
1. Introduction 

Ultracold quantum gases - Bose-Einstein condensates (BECs) - provide a uniquely 
accessible testbed for an astonishing range of effects occuring in diverse physical systems 
ranging from antiferromagnetic materials [1] and nonlinear photonic structures [2] to 
astrophysical black holes [3, 4]. The recent observation of a quantum three-body bound 
state in an ultracold cesium gas [5] is only one example of a fundamental effect that 
was predicted decades ago, but had not lent itself to experimental observation in any 
other physical system. The detailed exploration of these distinct physical phenomena 
in condensates is enabled by well-developed experimental techniques that provide 
unprecedented control over the interatomic interactions and the potentials confining the 
atoms. Advanced trapping schemes have recently been demonstrated by constructing 
box-like traps [6] and circular waveguides [7]. 

In this article we study a Bose-Einstein condensate with repulsive interatomic 
interactions confined in a harmonic trap and additionally exposed to a laser beam 
containing an optical vortex [8]. By tightly focussing the laser beam blue-detuned 
from an atomic resonance, it is possible to create a finite length optical tunnel for the 
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BEC, defined by the "dark" region of the beam in the vortex core. This dark region 
is located between the central vertical white lines in Fig. 1(a) and corresponds to the 
potential valley between the potential peaks in Fig. 2(a). 

The main purpose of this paper is to show that the proposed trapping configuration 
can enable the controlled study of a complex topological structure consisting of multiple 
stationary vortex rings in a BEC. Unsupported ring vortices are energetically unstable 
against contraction [9] and have so far been only observed in non- stationary situations, 
for example as the decay product of dark solitons or in soliton collisions [10, 11]. In 
contrast, our method could enable the stabilisation of a stack of up to nine parallel 
vortex rings under realistic experimental conditions. The creation of stationary vortex 
rings would also aid studies of their interactions [12] and represent a stepping-stone 
towards the creation of even more complex matter-wave structures such as skyrmions 
[13, 14, 15, 16]. 

A skyrmion is a topologial soliton in a two component BEC, in which the contraction 
of the ring vortex in one hyperfine component is prevented by filling it with a second 
hyperfine state carrying a line vortex. This also represents a tool to stabilise multiple 
ring vortices [26]. The optical tunnel system presented here results from a skyrmion if 
one replaces the hyperfine component carrying a line vortex by an optical field. Our 
system can thus be thought of as "atom-light-skyrmion" . 

Furthermore, in the proposed trapping configuration the BEC is flowing through 
the optical tunnel with velocities close to the speed of sound. Therefore we discuss 
the prospects of using our system for studies of analogue gravity effects. It has been 
discovered by W. Unruh [3] that within the hydrodynamic approximation [17], the 
equations of motion for quantised sound waves in a flowing BEC are analogous to 
those of a scalar quantum field in curved space time. Within this analogy, a sonic 
horizon, defined as a surface on which the normal component of the fluid flow becomes 
supersonic, corresponds to the event horizon of a black hole. Unruh proposed that 
a sonic horizon created in a BEC flow could enable the experimental observation of 
Hawking radiation [18, 19] of sound waves in a BEC. 

Here, we show that the BEC flow in the optical tunnel can be driven into the 
supersonic regime by inserting a laser generated "hump" potential into the tunnel. The 
perturbation can be engineered to result in the development of a single grey soliton in 
the condensate flow, which is immobile in the lab-frame. The BEC flow past the soliton 
is supersonic. 

The rest of this paper is organised as follows: In section 2.1 we introduce our 
model and describe the optical vortex potential. In section 2.2 we explain how the 
total potential has to be adjusted to allow stabilised ring flows, while we present the 
initial state of our imaginary time method in section 2.3. In section 3.1 we show some 
exemplary stable ring flows found numerically, followed by a brief description of a scheme 
for their experimental realisation in section 3.2. Section 3.3 contains a discussion of the 
dynamics that arise, when the flow is driven into the supersonic regime. Finally, we 
devote section 3.4 to an analysis of supersonic flow past a grey soliton in the context of 
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Figure 1. Cross-sections of density and velocity structures for different stationary 
ring flows, (a) Shows the BEC density (red highest, blue lowest) in the rz plane for 
case (i), defined by parameters: Vo = 4SHu) r , wq = 3/im, q = 7. The arrows indicate 
the BEC flow, thin white lines are isocontours of the optical vortex potential. The low 
isocontour (marked by L) corresponds to V v = 1 x 10~ 32 J, the high one (marked H) 
to V v = 1 X 10~ 31 J. The intersection of ring singularities with the rz plane is marked 
by x (o) for mathematically positive (negative) circulation. 

(b) BEC density and (c) Mach number M at r — along the z axis, (d) BEC 
amplitude in radial direction at z = 0. (e) BEC amplitude in radial direction at 
z = 170/im. Results for case (i) are indicated by solid lines in (b-e). Results for case 
(ii): Vq = 99hiu r , wq — 2.8/im, q = 9 are indicated by dashed lines in (b-e). 



analogue gravity. 



2. Model 



2.1. Gross-Pitaevskii equation 

We describe the dynamics of the BEC wavefunction ip with the Gross-Pitaevskii equation 
(GPE) [20]: 



dip 



(i) 



Using cylindrical coordinates r, z, ip, we model a BEC cloud in a cigar-shaped trap 
with frequencies u r radially and uo z longitudinally. Throughout the paper we assume 
that ip does not depend on (p, hence A = ^ + + The interaction coefficient 
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U is related to the s-wave scattering length a s : U = 4nh 2 a s /m, where m is the atomic 
mass. Both parameters are taken for 87 Rb: m = 1.44 x 10~ 25 kg and a s = 5.67nm. 
If the wavefunction is written in the polar form if) = y/pexp (m?), we can extract the 
atom number density p and velocity v = KS/'d/m. For later reference we also define the 
flow speed v = |v|, the local speed of sound c = ^JUp/m, the Mach number M = v/c 
and the healing length £ = h/(y/2mc). Finally V = V t + V v + V h denotes the external 



potential whose individual components are: 

V t (r,z) = ^m(cu 2 r 2 + u 2 z 2 ), (2) 

V v (r,z) = Vy l a (z)-^w 2 e W (-r 2 / a (z) 2 ) , (3) 

V h (r,z) = V ht cr t (z)- 2 w 2 exp (-r 2 /a t (z) 2 ) , (4) 

i=l,2 

a^z) 2 =w 2 + (z-zfymwi). (5) 



The repulsive potential V v can be generated by a blue-detuned optical vortex laser 
beam [21]. The charge of the optical vortex / is taken to be equal 2. In the dynamical 
simulations we employ a "hump" potential Vh, generated by one or more additional laser 
beams without any vorticity. They will be used to create small scale modifications of 
the overall potential experienced by the atoms within the core of the optical vortex. The 
propagation direction of the laser beams coincides with the long trap axis z. Further, 
wo,i,2 are parameters describing the waist of the laser beams at their focus, /co,i,2 are the 
wave numbers of the light and the zf* control the position of the beam focus. For the 
results presented here we used fcj = 1.4 x 10 7 m _1 and zf^ = 0. 



2.2. Potential design and trial states 

Figure 2 shows a slice in the rz plane of the combined potential V of the harmonic trap 
and optical vortex V = V t + V v . An equal energy contour is indicated by the white 
lines. In three dimensions the tight focus of the laser beam generates a circular local 
maximum, whose intersect with the rz plane yields the peaks situated symmetrically 
around the origin at \r\ ~ 5/im, see figure 2. For a suitable choice of parameters the 
combined potential exhibits saddle points, marked S in figure 2(a). If the value of 
the chemical potential is above the potential value at the saddle points, the BEC can 
classically "wrap" around the intense region of the optical vortex. 

We find parameter sets for which the saddle points exist using the following method: 
First we determine the position of the maximum of the vortex potential V v Eq. (3) with 
respect to r, keeping z fixed. We thus solve 

£v v (r,z) = (6) 
for r. We ignore the harmonic trap in this step. The solution of Eq. (6) is: 



= + ^ (7) 
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Figure 2. (a) Combined potential V of the harmonic trap Vt and optical vortex laser 
beam V v for Vo — 99hu> r , wq = 2.8/im, q = 9. Trap frequencies: to r — 7.8 x 2ir rad/s 
and u z = 0.5 x 2tt rad/s. Saddle points of the potential are indicated by S. The 
white line is an isocontour V = fiQ at p,o/(huj r ) = 20.9. (b, c) Radial cross sections 
of the potential at z — (b) and at z = — 210/im (c). Dotted lines indicate /io- (d) 
Dependence of the height of the potential "ridge" (local maximum) of 7 on z. This 
corresponds to the function V(r m (z), z) as described in the text. For \z\ > 280/im the 
local maximum disappears for the parameters chosen. 



Now taking Vt into account, the height of the potential "ridge" visible in figure 2(a) 
is approximately given by V(r m (z), z), where V t is now taken into account. Figure 2 
(d) shows V(r m (z), z). If saddle points exist, they must be minima of this function. 
For V(r m (z),z) to have extrema besides the obvious maximum at z — 0, the quartic 
expression 



d Tr/ . , A . 2e H 1 V Wq ( 4 z 2 \ IrrvJl 

must have real roots. After selecting the optical vortex parameters and u r , equation (8) 
yields a useful approximation for the strongest axial confinement uj z allowed. 



2.3. Energetically stable flows 

In order to find energetically stable ring flows we solve Eq. (1) for V = Vt + V v in 
imaginary time [22] and obtain stationary states of the system. The initial trial function 
for the imaginary time algorithm must have the same topological structure as the target 
ring flow. To this end, we multiply an amplitude a/ p(r, z) obtained from the Thomas- 
Fermi approximation [17] with a phase factor exp (ii](r, z))\ 

ip(r, z, 4>) = y/p(r,z) exp (ir](r, z)). (9) 

The phase is given by 

/ \ / \z\ (D - r)\ 
r](r, z) = q arccos I — . (10) 



z s(r, z) 
The quantity 

s( r ,z) = ^z 2 + (D-r) 2 (11 
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measures the distance to the single circular ring singularity implicit in Eq. (10). On 
any closed loop threaded through this singularity, the phase 77 varies by 27rg, where q 
is the ring vortex charge. The singularity is located at the point r = D, z = in our 
coordinates. An energetically stable excited state with high-charge ring flow corresponds 
to a local minimum of energy. If the trial condensate function is chosen to have the 
same topology as the energy minimum, we obtain the stabilised flow after a sufficiently 
long imaginary time [22]. 

3. Results 

3.1. Stationary states 

The low density region of a vortex ring in a BEC has a toroidal shape, which matches 
the shape of the high intensity region of a focussed optical vortex laser beam. In the 
presence of the optical vortex potential it becomes energetically favourable for the BEC 
ring singularity to be located at the maximum of laser light intensity. This prevents the 
contraction of the vortex ring, a mechanism similar to the stabilisation of line vortices 
within a harmonically trapped BEC cloud [25]. 

Examples of stationary states resulting from the imaginary time evolution are shown 
in figure 1. In a harmonic trap with uj r = 7.8 x 27rHz and oj z = 0.5 x 27rHz the BEC 
contains 2.2 x 10 6 atoms. The subplots are for two different sets of optical vortex 
parameters and ring flow charges: case (i) and case (ii), as given in the caption. Case (ii) 
represents the more intense optical vortex. It can be seen in figure 1(a) that the single 
charge 7 ring singularity that was present in the initial state has broken up into a regular 
stack of seven unit charge ring singularities. For case (i) they are responsible for the 
multi-peak structure of the Mach number shown in figure 1(c). If the parameters are 
altered towards those of case (ii), the influence of the singularities on the flow in the 
tunnel is reduced. This corresponds to the transition from the solid to the dashed lines 
in figure l(b-c). 

The integral along the z axis of the flow velocity is directly connected to the ring 
vortex charge q. Due to the quantisation of circulation in a BEC, its velocity must 
obey § c v ■ dl = ^q, where C denotes any closed contour threading through all the 
ring singularities. As the BEC has a much smaller cross-sectional area available when 
it flows through the tunnel than when it returns on the outer shell of the cloud, the 
velocity in the tunnel is much higher. The above equation thus approximately becomes: 



where L denotes the length of the optical tunnel in the z direction. 

The presented stability analysis of stationary states is not completely general, as we 
impose cylindrical symmetry throughout the imaginary time evolution. An important 
class of ring vortex excitations, the Kelvin modes (see Refs. [23, 24] and references 
therein), could in principle break this cylindrical symmetry. However, a full three- 




L/2 



(12) 
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dimensional (3D) stability analysis of skyrmions [15] indicates that the 2D model should 
give a very good indication of the stability properties: In Ref. [15], we have not found 
an indication of energetically unstable Kelvin modes of skyrmion's ring singularities. 
Nonetheless we point out that 3D studies of the optical tunnel would be an interesting, 
numerically challenging extension of our present work. 

3.2. Experimental creation scheme 

The method of stabilising the ring singularities by pinning, proposed in the present 
article, should be experimentally realisable with present technology. The procedure to 
form the ring flow consists of initially creating a BEC at rest, already in the presence of 
the optical vortex. Then a ring flow structure with high vorticity is seeded by means of 
phase imprinting. In Ref. [13] J. Ruostekoski and J. R. Anglin show in detail how the 
phase structure of the vortex ring, similar to the one given by Eq. (10), can be created. 
Their method uses a superposition of coherent light fields. In the presence of dissipation 
the BEC will, after the seeding, evolve towards stationary states like those in figure 1. 

Compared to the previously known method to stabilise the ring vortex as a part of 
a skyrmion, our proposal is more flexible. A skyrmion is a topological soliton in a two 
component BEC, composed of ring vortex in one hyperfine component, which is filled 
with a second hyperfine component carrying a line vortex. This filling prevents the 
contraction of the ring vortex in the skyrmion. [14, 15]. However, even the stabilisation 
of a charge two ring flow in a skyrmion requires large atom numbers of the order 1 x 10 7 , 
and more for higher charges [26]. In comparison, case (ii) represents a charge 9 ring flow 
in a condensate of only 2.2 x 10 6 atoms. 

3.3. Supersonic flow 

Supersonic flow is expected to be energetically unstable [27]. This manifests itself 
in the imaginary time evolution as a decay of the total ring flow charge from q\ to 
<?2 < Qi, whenever the parameters are such that a final state with charge q\ would 
breach the speed of sound. Despite the energetic instability, it is not generally true that 
the supersonic flow is also dynamically unstable [28] . Previously, it has been shown that 
dynamically stable flows with a supersonic region exist in a ID ring system [29, 30], and 
that they can be obtained from the subsonic ones by an adiabatic change of trapping 
potential parameters. 

Here we use the following sequence for the creation of supersonic flows: Initially the 
BEC is prepared in the stable state with subsonic ring flow corresponding to case (ii). 
Then an additional potential as in Eq. (4) is added dynamically. We found that for all 
potential widths, locations and ramp-up times employed, the generic response of the 
BEC was dynamical instability and the emission of a single or multiple grey solitons, 
travelling in the narrow confines of the optical tunnel. Refs. [31, 32] report on similar 
results in ID simulations. For our 2D simulations (that pertain to a 3D situation) 



Supersonic optical tunnels for Bose-Einstein condensates 



8 



we have employed an adaptive Runge-Kutta-Fehlberg method within the high-level 
programming language XMDS [33]. 

According to hydrodynamic theory, a hump potential with a single peak does not 
allow a quasi-stationary flow that is subsonic in a finite interval and subsonic on either 
end of this interval [34, 35]. Such a flow requires the realisation of a double de-Laval 
nozzle, using two hump potentials. The first of these renders the flow supersonic, while 
the second decelerates it back to subsonic velocities. We have numerically examined a 
variety of dynamical sequences, inserting double de-Laval potentials into the existing 
ring flows, and found that the result was always the same: The speed of sound is 
breached initially at the hump farthest in the direction of BEC flow, thus not realising 
a double de-Laval structure. Nonetheless, we can not exclude that with more careful 
adjustment of the potential the stationary ring flows presented here could be used as a 
staging point for double de-Laval nozzles. These have been suggested for the creation 
of sonic horizons in the context of analogue gravity [35, 36]. 

We now present our results achieved with the dynamical addition of a particularly 
shaped single peak potential. For the sudden addition of a hump potential, shaped as in 
figure 3(d), we find the following: Once the flow exceeds the speed of sound a single grey 
soliton is emitted, and for time spans as long as 0.58s it remains located near z — 0, see 
figure 3. Only after that time does the grey soliton begin to accelerate. The potential 
was obtained using Eq. (4) with parameters Vhi = —8.5hw r , Vh2 = 12.5ftw r , w\ = 1.2/xm 
and w 2 = 0.7/zm. 

We observe that at all times the flow speed in the tunnel is a function of z only, 
i.e. it remains quasi one dimensional. It is thus justified to apply ID soliton theory. 
We find excellent agreement between the functional form for ID grey solitons [37, 38] 
and the structures generated dynamically in our simulations. Further the dynamical 
behaviour of the soliton after its creation can be well understood with the aid of a ID 
Lagrangian approach [37, 38]. To this end, we apply the Ansatz [37, 38] 



Here fj, is the chemical potential, u(z,t) the soliton wavefunction and Pb s (z) denotes the 
background condensate density. We work in a frame of reference that moves with the 
speed of the background flow (i>b g ) at the soliton centre (z ). This speed is allowed to 
vary slowly in time. The parameters in the wavefunction A(t) and B{t) are related to 
the solitons width and depth and denotes the codensates healing length at z (t). 
Further we have A(t) + B{t) = 1. The leading order result for the soliton equation of 
motion following [37, 38] is: 



We find that this equation describes the trajectories of solitons in our simulations well, 
within uncertainties arising from the numerical extraction of pbg and v b g . In particular, 
Eq. (15) explains the immobility of the soliton shown in figure 3(a-b). In the vicinity 




(13) 
(14) 



u(z,t) = B(t) tanh [B(t)(z - z (t))/£(t)} +iA{t). 




(15) 




Figure 3. Scenario of long-time immobile soliton as sonic horizon, following a sudden 
change in the potential as shown in (d). (a) Time evolution of the Mach number M. 
Subsequent samples shifted by 2. Dashed horizontal lines indicate M = 1 for each 
time sample. The peaks of the Mach number are truncated at M = 2. (b) Position 
of the maximum of M. (c) (red) Shape of the ergoregion [39] and (black) the BEC 
cloud (shown is an isocontour at p = 3 x 10 16 m~ 3 ). (d) The potential is suddenly 
modified (at t = 0) from the shape indicated by the dashed line to the solid line. 



of the soliton's initial position the effective potential for the solitons, given by the 
background density, is approximately flat. Equation (15) suggests that the final travel 
direction of the soliton can be manipulated between positive and negative z by using 
potentials which are slightly asymmetric along that axis, which we have confirmed 
numerically The potential shape in such a case is still similar to figure 3 (d), but 
the wells adjacent to the local maximum are of different depths. The BEC density 
increases on the side with the deeper well, making it the favoured escape direction for 
the grey soliton. 

3.4- Prospects of analogue gravity 

The purpose of the present section is to consider whether and how grey solitons could 
play a role in studies of analogue gravity. Our goal here is not to give a rigorous and 
complete account, but to present basic ideas and stimulate further interest. 
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The grey solitons obtained in section 3.3 possess a sonic horizon where the BEC flow 
past them becomes supersonic, termed black hole horizon (BH). They also carry a white 
hole horizon (WH), the location where the flow returns to subsonic velocities once it 
has passed the soliton. The separation d between the black-hole and while-hole horizons 
depends on the soliton size, which is determined by the healing length. We numerically 
find that d ~ 2£. The conventional derivation of analogue Hawking radiation [4] requires 
the length scale of variations in the BEC flow, and hence the separation between BH 
and WH, to be much greater than the condensate healing length. The reason for this is 
that only modes with phononic wave numbers k such that k^ ^> 1 propagate in direct 
analogy to scalar fields in curved space-time [3, 4, 29, 40]. $. 

There exist interesting phenomena in analogue gravity, which crucially involve 
white holes in addition to black holes (e.g. the black hole laser effect [43]). However, 
BEC flow configurations with a widely separated black- and white-hole pair are usually 
dynamically unstable [29, 44], resulting in grey soliton emission from the white- hole [30]. 
This is also the origin of the solitons discussed in section 3.3. As a grey soliton already 
represents the end-product of dynamical instability, the black and white holes contained 
within it are stable, for as long as the soliton exists. Given the instabilities associated 
with well separated BH-WH pairs, compared to the ease with which one can obtain grey 
solitons, it would be interesting to enquire, which features of the Hawking or black-hole 
laser effects persist in a rapidly varying flow. Further motivation for research on quantum 
emission properties of rapidly varying condensate flows arises from the strict limits on 
the achievable analogue Hawking temperatures as long as the flow varies slow enough 
to be considered in the hydrodynamic regime [45]. It was shown that these arise due to 
three-body losses. 

It has been stated that the core ingredients of Hawking radiation are exclusively the 
existence of an apparent horizon and the validity of the eikonal approximation there [46]. 
The Mach number around the long time immobile soliton becomes as high as 600§. The 
group velocity v g of Bogoliubov excitations with wave number k is given by [40] : 



Inserting the speed of sound near the soliton, we see that only wave packets 
comprised of wavelengths shorter than g£ with g = 0.015 have a group velocity in 
excess of 600c. There is thus a window of wave numbers, 7r/£ <C k <C 2tt/ (<?£), for which 
a point of no return exists and the eikonal approximation should be valid in its vicinity. 
The location of this "horizon" varies depending on the major wavelength of the wave 
packet. Due to the steep Mach number profile, the "horizons" for different wavelengths 

| This issue is however a subtle one: All excitations experience a formally infinite blucshift i.e. increase 
in wave number, when traced back in time towards the horizon. They thus necessarily enter the non- 
phononic part of the Bogoliubov spectrum. However, it was argued in Rcfs. [41, 42] that the analogue 
Hawking effect persists despite this fact. 
§ The Mach number is truncated at M = 2 in Fig. 3(a). 



c 2 + ^k 2 



k 



(16) 
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are spatially close. Quantum emission of quasi-particles thus appears conceivable for 
wave numbers within the window, however further research is required to determine the 
existence of an effect. 

Another interesting analogue gravity system, for which a BH-WH pair associated 
with the propagation of a soliton was previously discussed is fermionic superfluid 3 He- 
A [47]. There the soliton of choice is a domain wall between regions of parallel and 
anti-parallel orientation of the spin and orbital angular momentum of Cooper-pairs. 
In such a system the separation between BH and WH can significantly exceed the 
coherence length [48]. Owing to the structure of the superfluid order parameter near 
the domain wall, the scenario discussed in Ref. [47] realises a rotating analogue BH, 
with an ergoregion boundary distinct from the apparent horizon [47] , in contrast to the 
simpler setup under consideration here. 

Returning to our condensate scenario, the velocity of the bulk BEC beyond the 
tunnel is effectively zero, see Fig. 1(c). This provides a region of flat space for the 
detection of quasi-particles emitted by the horizon. The steering of the soliton's escape 
direction could be an advantage for spectroscopic detection of horizon radiation [49]. 
The horizon-soliton can be made to travel towards z > 0, while the generated quantum 
excitations might be spectroscopically detected in the bulk BEC at z < before the 
soliton reaches the cloud edge. 

One of the advantages of the analogue gravity program over the astrophysical 
original, is the ability to push the studied system to the validity limits of existing 
derivations of horizon radiation. This allows the segregation of essential premises from 
merely analytically convenient assumptions. In this sense our system might one day 
complement more conventional analogue gravity setups along the lines proposed in 
Refs. [29, 50, 35, 36, 30]. 

4. Conclusions 

We have introduced a novel arrangement for the studies of persistent ring flows in a 
BEC. Our method utilises optical vortices, laser beams with a phase singularity at the 
centre and hence a zero intensity hollow core. If tightly focussed, the highest intensity 
region of the optical vortex can fit into the low density region of the matter-wave ring 
vortex in the BEC, stabilising the ring singularity against contraction. We numerically 
showed that energetically stable BEC states with ring flows of charges up to q = 9 exist 
in our system. 

We have further analysed the response of the ring flows to dynamically varied 
potentials which drive the flow into the supersonic regime, realising a supersonic optical 
tunnel for the BEC. In this case, our simulations show the creation of quasi one 
dimensional grey solitons in the tunnel. The solitons appear due to the dynamical 
instability once the speed of sound is reached by the velocity of the BEC flow. A single 
grey soliton that is immobile in the lab-frame for a long time was found for a particular 
choice of potential. 
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Finally, we have discussed the prospects of using this single, stabilised soliton as a 
tool for analogue gravity. We point out that whilst the Mach number around the realised 
horizon varies too rapidly to expect the usual form of analogue Hawking radiation, the 
soliton's sonic horizons might be suited to probe the physical boundaries for quantum 
emission from a horizon. 
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